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Observability of state variables and parameters of a dynamical system from an observed time 
series is analyzed and quantified by means of the Jacobian matrix of the delay coordinates map. 
For each state variable and each parameter to be estimated a measure of uncertainty is introduced 
depending on the current state and parameter values, which allows us to identify regions in state 
and parameter space where the specific unknown quantity can (not) be estimated from a given time 
series. The method is demonstrated using the Ikeda map and the Hindmarsh-Rose model. 


In physics and other fields of science including quanti¬ 
tative biology, life sciences, and climatology, mathemati¬ 
cal models play a crucial role for understanding and pre¬ 
dicting dynamical processes. In the following we assume 
that such a model exists and is known. But even in the 
ideal case of a model obtained from fundamental phys¬ 
ical laws this model typically contains some parameters 
whose values have to be determined depending on the 
physical context. Furthermore, not all state variables of 
the model may be easily experimentally accessible. To es¬ 
timate the unknown parameters and state variables you 
may either devise specific experiments focusing on the 
quantity of interest or you can try to extract the required 
information from a measured time series of the process 
to be modeled. Technically, several estimation methods 
exist, including observer or synchronization schemes 
E], particle filters [7], a path integral formalism , or 
optimization based algorithms [TOHIl]. However, these 
methods may fail and at this point the question arises 
whether the failure is due to the specific algorithm used 
or due to a lack of information in the available time series. 
In this article we address the second option and present 
a general approach for answering the question whether 
a given time series enables the estimation of parameters 
or variables of interest in a given model. The mathe¬ 
matical tool that is used to answer this question is delay 
reconstruction p^SHlT] and the basic criterion for local 
observability is the rank of the Jacobian matrix of the 
delay coordinates map. This approach was motivated by 
work of Letellier, Aguirre, and Maquet p^H2Q] who stud¬ 
ied the question which state variables can be estimated or 
observed from a given time series using derivative coordi¬ 
nates. Observability of (continuous) dynamical system is 
also a major issue in control theory [2Tti23] and nonlinear 
time series analysis [24]. Here we consider discrete time 
and delay coordinates, and we introduce a quantitative 
measure of uncertainty which in general varies on the at¬ 
tractor and thus indicates where in state space estimation 
is more efficient and less error prone. Furthermore, we 
focus not only on state variables but also on observability 
of model parameters. 

Let’s assume, first, that our model of interest is a M- 


dimensional discrete dynamical system 

x(n + l) =g[x(n),p] (1) 

given by an iterated function g depending on the state 
vector x(n) = (xi(n),... ,xm(^)) ^ at time n and 
K parameters p = (pi,... ,Pk) ^ • This system gen¬ 

erates the times series {s(n)} with s{n) = /i(x(n)) (for 
n = 1,..., A^), where h denotes a measurement or obser¬ 
vation function. The time series {s(n)} can be used to 
construct a D dimensional delay reconstruction HMT], 

y = (s(n), s{n + 1),...., s{n P D - 1)) (2) 

= G(x, p) G 

providing the delay coordinates map G : ^ 

To uniquely recover the full state x and the parame¬ 
ters p from the observations represented by the recon¬ 
structed state y the map G has to be locally invertible. 
More precisely, let M P K < D and let (x, p) eU where 
U C is a smooth manifold. Then G is locally in¬ 

vertible on the image G{U) C if the D x (M + K) 
Jacobian matrix UG(x, p) has full rank M -\- K (i.e., G 
is an immersion US]). 

The map from delay reconstruction space to the 
state and parameter space is locally given by the 

(pseudo) inverse of the Jacobian matrix DG of the de¬ 
lay coordinates map G, which can be computed using a 
singular value decomposition 

DG = USV^^ (3) 

where S = diag(cri,..., ctm+k) is a (M + K) x (M + K) 
diagonal matrix containing the singular values cri > <72 > 
... > cfm+k > 0 and U = ..., and V = 

(v^^\ ..., are orthogonal matrices, represented 

by the column vectors G and G re¬ 

spectively. is the transposed of V coinciding with 
the inverse V~^ = Analogously, = U~^ and 

the (pseudo) inverse Jacobian matrix reads DG~^ = 
VS-^U^^ where S-^ = diag(l/ai,..., l/aM+x). Mul¬ 
tiplying by U from the right we obtain DG~^U = VS~^ 
or 

{j = + K). (4) 

^3 
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FIG. 1. (Color online) The (pseudo) inverse Jacobian matrix 
DG~^{y) maps perturbations of y in delay reconstruction 
space to deviations from the state x whose magnitudes depend 
on the direction of the perturbation as described by Eq. 


In Fig. the transformation of singular vectors Eq. 
is illustrated for the case M = 2 and iF = 0 (no 
unknown parameters). The diagram shows how small 
perturbations of y in delay reconstruction space result in 
deviations from x in the original state space. Most rel¬ 
evant for the local observability of the (original) state x 
is the length of the longest principal axis of the ellipsoid 
given by the inverse of the smallest singular value <J 2 (see 
Fig. 0 . Small singular values correspond to directions in 
state space, where it is difficult (or even impossible) to 
locate the true state x given a finite precision of the re¬ 
constructed state y. The ratio cTmin/c^max of the smallest 
and the largest singular value is a measure of observabil¬ 
ity at the reference state x. By averaging on the attractor 
we define (analogously to a similar definition for deriva¬ 
tive coordinates mm) the observability index 


1 ^ 


- 

mm 

N ^ a 

n=l rnax 


(x)' 


(5) 


If the perturbations of y are due to normally dis¬ 
tributed measurement noise than they can be described 
by a symmetric Gaussian distribution centered at y 


Q{y) = 


exp[-|(y-yr-^S-Hy-y)] 

{2tt)^ det{T,y) 


(6) 


where y is the perturbed state, Sj, = diag(p^,..., p^) = 
P^Id denotes the Dx D covariance matrix {Id stands for 
the U-dimensional unit matrix), and the standard devi¬ 
ation p quantifies the noise amplitude. For (infinitesi¬ 
mally) small perturbations Ay = y — y this distribution 
is mapped by the pseudo inverse of the linearized delay 
coordinates map to the (non-symmetrical) distribution 


P(x) = 


exp 


-i(x-x)‘-'i]y(x-x) 


(7) 


v/(27r)^+^det(S,) 

centered at x with the inverse covariance matrix 
S-i = DG^’^T.pDG 
= jjDG^'^DG = 

The marginal distribution Pj of the jth state variable 
centered a Xj is given by 


(8) 


-g ) 


1 


pj\/^ 


exp 


(Xj 




2pj 


(9) 




FIG. 2. (Golor online) Observability of the state variables 
xi and X 2 of the Ikeda map Eq. 0 from 3i XI time series 
(with known parameters, iF = 0, M = 2). (a), (b) Golor- 
coded ratio of singular values cTmin/cTmax vs. xi and X 2 for 
reconstruction dimension D — 2 {b) and D — ?> (b). The 
white curves in (a) indicate the location of zeros of del{DG). 
(c), (d) Golor-coded uncertainties vi (c) and 1^2 (d) of xi and 
X 2 estimates, respectively. Note the logarithmic color axes. 
Black dots represent the Ikeda attractor. 


where the standard deviation pj is given by the square 
root of the diagonal elements of the covariance matrix 
pj = y^Tixjj that can be obtained by inverting 11“^ 
[given in Eq. Q]. Since the noise level p of the observa¬ 
tions appears in Eq. (§ as a factor only we can, without 
loss of generality, choose p =1 and use 

Vy = sJ[DG^^DG]-^ = ( 10 ) 

as a measure of uncertainty when estimating Xj^ which 
can be interpreted as a noise amplification factor. The 
same reasoning holds for the unknown parameters p. 

To illustrate this quantification of observability we 
first consider the Ikeda map [25] z{n 1) = pi -\- 
P 2 ^(^) exp[ip 3 — ip 4 /(l + |2:(n)p)] with z{n) = Xi{n) + 
ix 2 {n) G C that can also be written as 

xi(n-hl) = Pi+p2[^i(^)cos6>n-X2(n)sin6>n] 

X2(n-hl) = P2[^i(^) sin6>^ + X2(n) cos6>^] ^ ' 

where = Ps — P 4 /[l-hx^(n)-hX 2 (n)]. For the standard 
parameters pi = 1, P 2 = 0.9, ps = 0.4, and p 4 = 6 this 
map generates the chaotic attractor shown in Fig. 

First, we consider a case where all parameters are 
known and only the variables xi and X 2 have to be es¬ 
timated from the observable s{7i) = xi{n) (i.e., M = 2 
and iF = 0). Figures [^a) and [^b) show (color-coded) 
the ratio of the smallest singular value amin = cfm and 
the largest singular value cTmax = cri of the Jacobian ma¬ 
trix DG{x.) of the delay coordinates map vs. xi and X 2 . 
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Reconstruction dimensions are D = 2 in Fig. |^a) and 
D = 3 in Fig. |^b), respectively. For D = 2^ the white 
curves indicate the zeros of the determinant of p) 

that are computed as contour lines. As can be seen parts 
of the Ikeda attractor cross these singularity manifolds 
or are close to regions in state space where the ratio 
<^min/<^max IS vcry closc to zcro, indicating an almost sin¬ 
gular Jacobian matrix DG. There, state estimation is 
not possible, a fact that reconfirms previous results indi¬ 
cating that reconstruction dimensions D > 2 are required 
for the Ikeda map [26]. For D = 3 the singularities dis¬ 
appear and only some regions with relatively low ratios 
o-min/o-max remain. 

Figures [^c) and id) show ui and U 2 versus Xi and 
X 2 , respectively. For both variables their uncertainties 
Uk vary and there are regions of low ui but relatively 
large z^ 2 - 

Figures |^a) and ib) show histograms of and U 2 
for different reconstruction dimensions D which were ob¬ 
tained from an orbit of length N = 1000000 on the Ikeda 
attractor. Due to the choice s{n) = xi{n) the uncer¬ 
tainty ui of XI is for all dimensions equal or less than 
one. For D = 2 the uncertainty 1^2 of X 2 reaches very 
high values > 10 ^ when the orbit passes those regions in 
state space where the Jacobian matrix DG is (almost) 
singular [see Fig. |^a)]. For reconstruction dimensions 
D = 3 the z/ 2 -histogram is bounded by 1^2 < 10^ in¬ 
dicating a significant improvement and for D = 4 the 
bound reduces to z ^2 < 10 , a value that doesn’t change 
anymore if the reconstruction dimension is increased fur¬ 
thermore. This feature is in very good agreement with 
previous results obtained when estimating Lyapunov ex¬ 
ponents from Ikeda time series [26| . 

To obtain the histograms shown in Fig. and in the 
following figures the model equations are used to gener¬ 
ate a trajectory which provides a representative sample 
and subset of the attractor (similar to numerical compu¬ 
tations of Lyapunov exponents). 

For the results shown in Figs. I^a) and [^b) only the 
state variables are estimated and all parameters are as¬ 
sumed to be known (M = 2, LC = 0). Figure shows 
also the uncertainties z/ 3 , z/ 4 , z/ 5 , and uq of the parame¬ 
ters pi, P 25 Ps, and p 4 for an estimation task where all 
variables (M = 2) and all parameters {K = 4) are un¬ 
known. For increasing reconstruction dimension D the 
distributions of all uncertainties converge with monoton- 
ically decreasing upper bounds (largest z^-values quanti¬ 
fying large uncertainty of estimates at specific locations 
on the attractor). 

Delay reconstruction can also be applied to observables 
s{t) = h[x(t)] from continuous dynamical systems, 

x = f(x,p), (12) 

using a suitable delay time r: 

y = {s{t), s{t + r), s{t + {D- l)r)) = G(x, p) e 

The Jacobian matrix DG{k^ p) of the delay coordi¬ 
nates map G can be computed by solving linearized 



FIG. 3. (Color online) Histograms (color-coded) of uncer¬ 
tainties m (a) and 1^2 (b) computed from Si xi time series of 
length N = 1000000 generated on the attractor of the Ikeda 
map Eq. 0 with reconstruction dimensions ranging from 
D = 2 to D = 7. The state variables xi and X 2 are esti¬ 
mated (M = 2) while all parameters are assumed to be known 
(X = 0). 



FIG. 4. (Color online) Histograms (color-coded) of uncer¬ 
tainties of state and parameter estimates of the Ikeda map 
Eq. (11) for reconstruction dimensions ranging from D = 6 
to D = 12. Distributions are computed from a time series 
of length N = 1000000 generated on the Ikeda attractor. All 
variables (M = 2) and all parameters {K = 4) are assumed 
to be unknown. 


equations providing the Jacobian matrices Da,(/)^(x, p) 
and ^0^(x, p) of the flow (j)^ generated by the system 
Eq. ( [T^ [27] . To demonstrate the application of the pro¬ 
posed uncertainty analysis to continuous time system we 
use the Hindmarsh-Rose (HR) neuron model [28] 

Xl = -x\ + Pix\ -h X2 - X3 

±2 = 1 -P2x‘l - X2 (13) 

^3 =P3 (^1 +P4(P5 -^ 3 ))- 

For parameter values pi = 3, p 2 = 5, p 3 = 0.004, p 4 = 
3.19, P 5 = 0.25 the HR model exhibits chaotic bursting 
of Xl and X 2 and slow variations of X 3 m- 

Figures [^a) andf^b) show the dependence of proba¬ 
bility distributions (S)lor-coded) of uncertainties z/ 2 , and 
z/ 3 , respectively, on the delay time r chosen for perform¬ 
ing the delay reconstruction. The reconstruction dimen¬ 
sion equals D = 7. With this example, all parameters 
are assumed to be known {K = 0) and the first state 
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FIG. 5. (Color online) Probability distributions (color-coded) 
of uncertainties 1^2 and when esti mat ing the state variables 
xi, X 2 ^ and X 3 of the HR-model Eq. ( p^ from a xi time series. 
In (a) and (b) the delay reconstruction dimension is fixed at 
D = 7 and the delay r is varied, (c), (d) show distributions 
for T = 0.1 and different reconstruction dimensions D. Cor¬ 
responding columns (histograms) of all four diagrams show 
results for the same window in time {D — 1)t used upon de¬ 
lay reconstruction, (e), (f) Observability index 7 ® (solid 
curve), cTmin (dotted curve), and cTmax (dashed cur^ vs. r 
and vs. D. 
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FIG. 6 . (Color online) Distributions of uncertainties uj vs. 
reconstruction dimension D obtained for the HR model (13) 
where all three state variables and all five parameters are 
estimated from a time series. The delay time r = 0.1 is 
fixed. 


variable is chosen as measured time series s{tn) = xi{tn) 
with tn = nr. Therefore, the estimation of xi is not 
much affected by the choice of the delay time and < 1 
(with ui ^ 1 most of the time, not shown here). As 
can be seen the centers of both distributions decrease 
monotonically with r indicating an improvement of the 
estimation accuracy for larger delay times. Figures [^c) 
and id) show histograms (color-coded) of uncertainties 
z/ 2 , and z /3 versus reconstruction dimension T) for r = 0 . 1 . 
Larger D provides lower uncertainties Uj and compared 
to Figs. I^a) and ib) very large Uj do not occur any¬ 
more. Note that corresponding columns of Figs. ia) 
and[^b) and Figs. ic) and [^d), respectively, are com¬ 
puted using delay coordinates covering the same win¬ 
dow in time ranging from r{D — 1) = 0.1 • 6 = 0.6 to 
r{D — l) = 30.1-6 = 1806-0.1 = 180.6. The more densely 
sampling (r = 0.1) underlying Figs. |^c) and |^d) pro¬ 
vides more information about the underlying dynamics 
and results in lower uncertainty values. Figures [^e) and 
l^f) show the observability index 7 Eq. ^ and mean 
values of the smallest and the largest singular values 
(Jmin and cTmax vcrsus r and respectively. While 7 
exhibits a clear peak, cTmin converges to an asymptotic 
value, and (Jmax increases monotonically, i.e., the lengths 
of the ellipsoid axes in Fig. decrease (l/cTmax) or con¬ 
verge (l/cTmin). 

If in addition to the three state variables xi, X 2 , and 
xs also the five parameters ,..., ps of the HR-model 
Eq. (13) are to be estimated from the xi time series 
then we have to cope with an estimation task with 


M + iF = 3 + 5 = 8 uncertainties whose distributions 
for r = 0.1 are shown in Fig. [^for delay reconstruction 
dimensions ranging from D = S to D = 2008. For in¬ 
creasing D the uncertainties z/i,..., z/e corresponding to 
xi,X 2 ,X 3 ,pi,P 2 ,P 3 decrease to values close to or below 
one. The uncertainties uj and us of parameters p 4 and 
P 5 , respectively, remain rather large (> 1000 ) even for 
high dimensional reconstructions. This feature indicates 
that it is very difficult to estimate both parameters to¬ 
gether. In fact, if p 4 (or p^) is known and only (or p^) 
has to be estimated (together with xi^X 2 ^xs^Pi^P 2 iP 3 ) 
then the uncertainty values of p^ (or ^ 4 ) are much smaller 
and lie in the range of the uncertainties of the other pa¬ 
rameters. Applying a state and parameter estimation al¬ 
gorithm [II1I29] we also encountered problems (in terms 
of large deviations from the true values) when trying to 
estimate both parameters p^ and ps together. These two 
parameters are to some degree redundant in the sense 
that different combinations yield (almost) the same xi 
time series and thus cannot be clearly distinguished us¬ 
ing a time series, only. 

The presented approach for quantifying uncertainties 
of model based state and parameter estimation from time 
series provides a general criterion whether and how reli¬ 
ably specific model variables and parameters can be es¬ 
timated from time series. This method is independent 
from any particular estimation method and it can be ex¬ 
tended in several ways, including unknown parameters in 
the measurement function and multivariate time series. 
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High uncertainty implies that the corresponding quan¬ 
tity of the model has small impact on the output and 
may thus be a candidate for reducing the formal model 
complexity by pruning. Furthermore, the information 
provided by the values of uncertainty can be exploited to 
improve state and parameter estimation methods. 
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